Vorticity structuring and Taylor-like velocity rolls triggered by gradient shear bands 
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We suggest a novel mechanism by which vorticity structuring and Taylor-like velocity rolls can 
form in complex fluids, triggered by the linear instability of one dimensional gradient shear banded 
flow. We support this with a numerical study of the diffusive Johnson-Segalman model. In the 
steady vorticity structured state, the thickness of the interface between the bands remains finite in 
the limit of zero stress diffusivity, presenting a possible challenge to the accepted theory of shear 
banding. 

PACS numbers: 47.50.+d, 47.20.-k, 36.20.-r. 
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I. INTRODUCTION 

Many complex fluids exhibit flow instabilities that re- 
sult in spatially heterogeneous, shear banded states. Ex- 
amples include wormlike [J, H, Q and onion [!, 0, 
surfactants; side-chain liquid crystalline polymers 
viral suspensions H H[; telechelic polymers [l(J; soft 
glasses fill] ; polymer solutions [13] ; and colloidal suspen- 
sions [13j |. In many cases, the instability is explained by a 
region of negative slope dT xy /d^ < in the constitutive 
relation T xy (7) between shear stress and shear rate for 
homogeneous flow, as shown in Fig [T^,. In this regime, 
homogeneous flow is unstable 14j with respect to the for- 
mation of bands of differing shear rates 71 and 72 , with 
layer-normals in the flow-gradient direction. Force bal- 
ance requires that the shear stress T xy is uniform across 
the gap, and therefore common to both the bands. Any 
change in the overall applied shear rate 7 causes a change 
in the relative volume fraction / of the bands according 
toalever rule7 = /7i + (l — /)72, while 71,72 andX^ re- 
main constant. In bulk rheometry, this leads to a plateau 
in the steady state flow curve at some stress T xy = T* 
(Fig. QJ,), the value of which is selected by accounting for 
spatial non-locality in the constitutive dynamics of the 
viscoelastic stress [15] . In what follows, we shall refer 
to the effect just described as gradient shear banding, or 
simply gradient banding. 

In some systems, flow induced heterogeneity has been 
reported in the vorticity direction. By analogy with 
the above discussion, this effect is often termed vorticity 
banding. It has been observed in onion surfactants [1, Ha ] 
and colloidal crystals [H| , accompanied by a steep stress 
"cliff" in the flow curve. Multiple turbid and clear vor- 
ticity bands also occur in some polymeric fl2| and mi- 
cellar [l7l [l8j solutions, accompanied by shear thicken- 
ing. In these cases, the bands not only alternate in space 
but also oscillate in time. In shear thinning viral suspen- 
sions, multiple stationary vorticity bands can arise in the 
regime of isotropic-nematic microphase separation 0, Q . 

In comparison with gradient banding, vorticity band- 
ing is poorly understood theoretically. To date, the main 
attempts to model it have invoked the analogy with gra- 
dient banding [lj|. As discussed above, gradient bands 



have different shear rates 71 , 72 , and coexist at a common 
shear stress T*. For vorticity bands, the moving rotor im- 
poses a shear rate that is common to each band. Pursuing 
the analogy with gradient banding, it is natural to invoke 
an underlying constitutive curve of the form in Fig.[Tp, in 
the case of shear thickening systems. This allows bands 
(A and B) of differing shear stresses to coexist, having 
layer normals in the vorticity direction. In steady state, 
one then expects a flow curve with a steep stress cliff, 
consistent with the experimental observations discussed 
above 0, [H, [l6| ■ This scenario can be adapted to shear 
thinning systems by invoking a constitutive curve of the 
form sketched in Fig. [IF,. In fact, at the level of ID calcu- 
lations performed separately in the flow-gradient and vor- 
ticity directions, such a curve can support either gradient 
or vorticity banded states [2(| • Which of these (if either) 
would be selected in a full 3D calculation remains an 
outstanding question. Indeed, any concrete calculation 
of vorticity banding to date has taken a simplified one- 
dimensional (ID) approach, allowing spatial variations 
only in the direction of the layer normals, and thereby 
imposing axial and radial symmetry. 

Beyond the traditional shear banding literature, other 
flow instabilities are well known to trigger structuring in 
the vorticity direction, in both simple and complex fluids. 
For simple liquids in curved Couette flow, the unstable 
centripetal interplay of fluid inertia with cell curvature 
gives rise to Taylor velocity rolls [2l[ stacked in the vor- 
ticity direction. An analogous inertia-free instability oc- 
curs in viscoelastic fluids, here triggered by viscoelastic 
hoop stresses [13, HH , and again leading to velocity rolls 
stacked along the axis of the Couette cylinder. The rolls 
are typically observed as bandlike structures, imaged by 
seeding the fluid with mica flakes. 

In contrast to the ID vorticity banding scenario dis- 
cussed above, both traditional and viscoelastic Taylor- 
Couette instabilities arc inherently 2D (at least), with 
velocity rolls comprising a circulation of fluid in the 
flow-gradient /vorticity plane [H|,[13]. Not unexpectedly, 
given this roll-like structure, the wavelength of the asso- 
ciated banding in the vorticity direction is roughly set by 
the width of the gap in the flow-gradient direction [2l| . 

In view of the above discussion, it is natural to ask 
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FIG. 1: a) Homogeneous constitutive curve and steady state 
flow curve for ID planar gradient banded flow in the DJS 
model at a = 0.3,7) = 0.05. Inset: schematic arrangement of 
the bands in curved Couette flow, b) Schematic constitutive 
curve giving shear thickening vorticity banding, c) A shear 
thinning curve allowing gradient or vorticity banding in ID. 

whether any link exists between traditional (ID) vorticity 
banding and 2D viscoelastic Taylor-Couette instabilities; 
or whether the two effects are entirely distinct. Given 
that both can be accompanied by shear thickening in the 
bulk flow curve, this is a difficult question to address ex- 
perimentally. Even if they are distinct in theory, it seems 
feasible that some experimental observations that have 
traditionally been interpreted as ID vorticity banding in 
fact comprise 2D viscoelastic Taylor Couette rolls. Likely 
candidates include those systems in which the wavelength 
of the alternating vorticity bands is comparable to the 
width of the cell in the flow-gradient direction, suggest- 
ing an underlying roll structure @, H El ES EH ■ This 
was recently suggested in the context of viral suspensions 
in Ref. 24]. In other systems, particularly those showing 
a very marked stress cliff in flow curve, the traditional 
ID scenario of Fig. [TJd or c remains more likely. 

In this paper, we suggest a novel mechanism by which 
vorticity structuring can emerge in complex fluids. A key 
feature of our approach is that, to some extent, it uni- 
fies traditional ID (gradient) banding descriptions with 
those of 2D roll-like instabilities. The basic idea is as fol- 
lows. The constitutive curve of Fig. QJi gives rise initially 
to ID gradient bands, via the conventional instability in 
the region of negative slope, dT xy jdrf < 0. These then 
undergo a secondary linear instability 0, |2f3], due to 
the action of normal stresses across the interface between 
the bands. This leads finally to pronounced undulations 
along the interface, with wavevector in the vorticity direc- 
tion. These are accompanied by 2D Taylor-like velocity 
rolls stacked in the vorticity direction, and undulatory 
vorticity stress structuring superposed on the underly- 
ing gradient bands. In contrast to the conventional in- 
ertial [2l| and viscoelastic [22j Taylor mechanisms, the 
vorticity instability introduced here does not rely on cell 
curvature, but occurs even in the limit of planar shear, 



to which our calculations are confined for simplicity. 

The results to be presented are in good agreement with 
recent experiments in which a gradient banded solution 
of wormlike micelles was found to be unstable with re- 
spect to interfacial undulations with wavevector in the 
vorticity direction [27j ■ We will return to a detailed com- 
parison with these experiments later in the manuscript. 
Our results might also apply to systems in which shear 
thickening [H and/or vorticity banding pj, EE E3, El 
is reported to set in at the right hand edge of a stress 
plateau in the flow curve (suggestive of underlying gra- 
dient banding, as in Fig. Hk); or in which gradient and 
vorticity banding have actually been observed concomi- 
tantly [291 ]. 

The paper is structured as follows. In Sec. |TT] we in- 
troduce the diffusive Johnson Segalman model, within 
which the subsequent calculations are to be performed. 
In Sec. IIIII we discuss ID calculations, confined to the 
flow-gradient direction. These predict gradient banding 
for applied shear rates in the regime of negative slope in 
the homogeneous constitutive curve. In Sec. lIVI we switch 
to two dimensions - the flow-gradient/ vorticity plane - 
and show an initially ID gradient banded "base state" 
to be linearly unstable with respect to undulations along 
the interface with wavevector in the vorticity direction. 
We then perform a full 2D nonlinear numerical study of 
the subsequent growth and eventual saturation of these 
undulations. Details of the numerical method are dis- 
cussed in Sec. [Vj followed by presentation of the results 
in Sec. lVIl Finally we summarise our findings and discuss 
some directions for future study. 



II. MODEL AND GEOMETRY 

The generalised Navier-Stokes equation for a viscoelas- 
tic material in a Newtonian solvent of viscosity r\ and 
density p is 

p(d t + v.V)v = V.(£ + 2r?D - PI), (1) 

where v(r, t) is the velocity field and D is the symmet- 
ric part of the velocity gradient tensor, (Vv) a( g = d a vp. 
The pressure field P(r, t) is determined by enforcing in- 
compressibility, 

V-v = 0. (2) 

The quantity S(r, t) in Eqn. Q] is the extra stress con- 
tributed by the viscoelastic component. In principle, the 
dynamics of this quantity should be explicitly derived by 
averaging over the underlying microscopic dynamics of 
the viscoelastic component. This was done in Refs. [30I ] 
for wormlike micelles. For simplicity, however, we use the 
phcnomenological Johnson-Segalman (JS) model [3l[ 
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(d t + v • V) £ = o(D • £ + £ • D) + (£ ■ ft + n ■ £) + 2GD + — V 2 £. (3) 



In this equation, G is a plateau modulus, r is the vis- 
coelastic relaxation time, and $1 is the antisymmetric 
part of the velocity gradient tensor. For a — 1 and 
£ = 0, Eqn. [3] reduces to the Oldroyd B model, which 
can be derived by considering the dynamics of an ensem- 
ble of Hookean dumbbells in solution. For \a\ < 1, the 
JS model captures non-affine slip between the dumbbells 
and the solvent, leading to the drastic shear thinning 
of Fig. 0Jl. Accordingly, a is called the slip parameter. 
The JS model is the simplest tensorial model to exhibit 
a regime of negative slope in the homogeneous constitu- 
tive curve, and so to predict a shear banding instability. 
As discussed further below, the diffusive term V 2 £ in 
Eqn. [3] is needed to correctly describe the ultimate shear 
banded flow 

Within this model, we study planar shear between par- 
allel plates at y — 0, L, with the top plate driven at ve- 
locity Vx~ At the plates we assume boundary conditions 
of dyYiap = V a, (3 for the viscoelastic stress, with no 
slip and no permeation for the fluid velocity. In the lin- 
ear stability analysis of Sec. IIV1 we consider small values 
of the Reynolds number Re = pi? jr\. In the nonlinear 
study of Sees. |V] and IVI1 we set Re = at the outset. 
Throughout we use units in which G = 1, r = 1 and 
L = 1. 



III. ID GRADIENT BANDS 

As noted above, to capture shear thinning the DJS 
model invokes a slip parameter a with \a\ < 1, giving 
non affine deformation of the viscoelastic component [3l| . 
The homogeneous constitutive curve T xy = £ XJ/ (7, a) +777 
is then capable of non-monotonicity, as in Fig. [T^,. For 
an imposed shear rate 7 = V/L in the region of decreas- 
ing stress, homogeneous flow is unstable with respect to 
fluctuations with wavevector in the flow-gradient direc- 
tion y 32]. A ID calculation then predicts separation 
into gradient bands of differing shear rates 71,72, with a 
flat interface in between. The diffusive term in Eqn. [3] is 
needed to account for spatial gradients of the shear rate 
and viscoelastic stress across the interface, which has a 
characteristic thickness 0(f). It also ensures a unique, 
history- independent banding stress T xy = T* [15| , as seen 
experimentally. We expect I = O(10 -4 ), set by the typi- 
cal micellar mesh size, in units of the (typical) gap size. 

IV. LINEAR INSTABILITY 

In Refs. [1^, (2(|, we considered the linear stability of 
this ID gradient banded state with respect to 3D (x, y, z) 
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perturbations of infinitesimal amplitude. In the flow di- 
rection x and vorticity direction z these are decomposed 
into Fourier modes with wavevectors q = g^x. + q z i. 
(Ref. [33[ had previously considered q = q x x in the 
pathological limit £ = 0, assuming "top jumping" .) For 
diffuse interfaces, I ~ 0.015, the ID state is linearly sta- 
ble. For I $ 0.015, we find it to be linearly unstable 
with respect to modes with wavevector q = q x yi in the 
flow direction (2f| [26[ . The associated eigenfunction es- 
sentially corresponds to undulations along the interface. 
For I ~ 0.005, the ID state is also unstable with respect 
to undulations with wavevector q = q z z in the vorticity 
direction. However, these modes are predicted to grow 
much more slowly than those with q = q x 5<L. Accordingly, 
in Refs. (25l . [26| . we focused mainly on the dominant 
modes, q = q x it. Subsequently in Ref. [34[, however, we 
showed that these undulations are cut off, once they at- 
tain finite amplitude, by the nonlinear effects of shear. In 
contrast, the vorticity direction is neutral with respect to 
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FIG. 2: a) Dispersion relation for perturbations about a ID 
banded state for £ = 0.00375, 0.00250, 0.00125. a = 0.3, 77 = 

0. 05, p/rj = 0.02. Inset: Linear dynamics of 2D code, starting 
from a flat interface. Solid lines: weight in modes, q z L z /2ir = 

1, 2; dashed: analytical prediction, b) Peak u* (q* in inset) in 
dispersion relation vs. 7 across the plateau of Fig.[]J,. Vertical 
dashed lines denote the edges of the stress plateau. 
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the shear. Accordingly, the modes with q = q z i should 
not suffer this cutoff and are therefore likely to contribute 
significantly to the ultimate nonlinear state, despite their 
much slower initial growth rate. With this motivation, 
in this paper we study the dynamics of the model in the 
flow-gradient /vorticity (y — z) plane. For simplicity and 
computational efficiency, we will assume uniformity in 
the flow direction x, returning in Sec. I VII below to com- 
ment on the validity of this simplification. It corresponds 
to taking a vertical slice through one side of an axisym- 
metric flow state in the planar limit of a Couette device. 

The growth rates u> of the vorticity modes q = q z z are 
shown in Fig. [2]}, at a single value of the imposed shear 
rate. States with thinner interfaces (smaller €) are more 
unstable (larger u> > 0). Fig.[2jD shows the growth rate of 
the maximally unstable mode for shear rates across the 
stress plateau of Fig. [1^. The corresponding wavelength 
A* = 0(1) is of order the rheometer gap L = 1 (inset). 
At small £, the instability persists across most of the 
plateau, so is likely to be unavoidable experimentally. 
The mechanism of instability is not fully understood, but 
is likely to stem from steep gradients in the normal stress 
and shear rate across the interface [25|, [3f| [36| . 

V. NUMERICAL METHOD 

To study the undulations once they have grown to at- 
tain a finite amplitude, beyond the regime of linear in- 
stability, we solve the model's full nonlinear dynamics 
numerically. In this section we discuss the details of 
our numerical method, which is adapted from that of 
Refs. [34|, |37]. Readers who are not interested in these 
issues can skip straight to Sec. IVII without loss of thread. 

The model equations have already been specified in 
Sec. HIl together with the flow geometry, boundary con- 
ditions and choice of adimensionalisation. For computa- 
tional efficiency, our numerical study is confined to the 
(y — z) plane, assuming translational invariance along the 
flow direction x. In the vorticity direction we take a cell 
of length L z , with periodic boundary conditions. 

We consider the limit of zero Reynolds number, in 
which Eqn. [1] reduces to 

= V.(S + 277D - PI). (4) 

To ensure that the incompressibility constraint of Eqn. [2] 
is satisfied always, we express the velocity in terms of 
stream functions <f> and ip: 

V x = dy(j), Vy = 8 Z tp, V Z = -Oylp. (5) 

In this way, Eqn. [2] need no longer be considered: it re- 
mains only solve Eqns. [3] and 01 with the velocity ex- 
pressed as in Eqn. [5] 



To solve these, the basic strategy is to step along a grid 
of time values t n = nAt for n = 1, 2, 3 ■ • ■ , at each step 
updating <f\ -0" -> S n+1 , n+1 , tp n+1 . Discretiza- 
tion with respect to time of any quantity / is denoted 
f(t n ) = f n , or sometimes below by f\ n . At each time- 
step, we first update the viscoelastic stress S™ — > 
using the constitutive equation |3|) with fixed, old val- 
ues of the stream- functions <j> n ,ip n . We then update 
(f> n ,ip n — > (f) n+1 , tp n+1 using the force balance equation 
(j4| with the new values of 

The update S™ — > S™ +1 using the viscoelastic consti- 
tutive equation ([3]) is performed as follows. As a prelim- 
inary step, we rewrite ((3J in the form 

3 ( E = /(Vv,S)-vVS + f¥S, (6) 

in which /(Vv, S) comprises the non-diffusive terms 
from the right hand side of Eqn. [3l In what follows, 
the three terms on the right hand side of Eqn. [5] are 
referred to as the local, advective and diffusive terms re- 
spectively. Numerically, they are dealt with in three suc- 
cessive partial updates S™ S™ +1 / 3 , S" +1 / 3 -> s n + 2 / 3 
and £™+ 2 / 3 -> S n+1 . 

In the first of these, the local term is handled using 
an explicit Euler algorithm [38|, checked for consistency 
against a fourth order Runge-Kutta algorithm 38] . Tem- 
porarily setting aside the issue of spatial discretization, 
the Euler algorithm can be written 

£ n+1/3 (?/, z) = S" + At /(Vv n , £"). (7) 

In terms of the stream-functions ip and 0, the velocity- 
gradient tensor Vv has Cartesian components 

/ \ 

vv = d 2 v <i> d y d z ip -d 2 v i> , (8) 

\ d y d z (j) d 2 z ip -d y d z i\) J 

in which we have omitted the superscripts n for clarity. 
Eqns. [7] and [8] are then spatially discretized on a rectan- 
gular grid in real space. In some runs of the code, the grid 
points are linearly spaced, with Zi = iAz for i = 1 ■ ■ • N z 
and yj = jAy for j — 1 ■ • ■ N y . In others, we used a 
nonlinear mapping in the y direction to focus attention 
on the region explored by the interface. For simplicity, 
most of the description of this section will concern the 
linear grid, though we will return briefly at the end of 
the section to discuss the nonlinear modification. In ei- 
ther case, any spatially discretized function / is denoted 
f(y J ,z l ) = fij, or sometimes f\ij. (The apparently un- 
usual order of the indices is a historical convention on 
the part of the author, stemming from a previous study 
in the x — y plane.) Eqn. [7] then becomes 
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S" +At/(Vv",£" 



(9) 



The derivatives in the components of Vv™ are discretized (in the case of a rectangular grid) as follows: 
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Ay 2 



(10) 
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4AxAy 



(ii) 



(12) 



Corresponding derivatives of <j) are obtained in the same way, replacing -0 by in the above equations. For values of 
ij at the edges of the flow domain, these formulae link to values of the flow variables at "phantom" grid points that 
lie just outside the domain. These values are specified by imposing the boundary conditions, the spatial discretization 
of which is discussed at the end of this section. 

The advective term is also handled using an explicit Euler algorithm (38|, on the same real space grid: 



■yn+2/3 _ y, 

ij ~~ V 



+ 1/3 
-1/3 



At kis^eis 



Atid^dy^-dy^d^) 

The derivatives of tp in this equation are discretized as follows: 

1 
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[ipi( j+ i) - &Cj-i)] and d *^\ij = [^(<+i)j - V>( 



-ibJ 



2Ay L^Wf^J r%v-x U ^-rnj 2Az 

The derivative of X with respect to y in Eqn. [13] was discretized using third-order upwinding 3£ 

1 



6Ay 



3sr 



2£ 



if w s |g->0, 



while 



(13) 



(14) 



(15) 



6Ay 



-i(j+2) 



6S 



i(j+l) 



3EJ 



2£ 



iO'-i) 



if wX-<0, 



(16) 



with analogous expressions for the derivative of S with respect to z. 

The diffusive term is handled by discretizing on y in real space as above, taking a fast Fourier transform z —> qi in 
the vorticity dimension using a standard NAG routine [40| , and solving the resulting problem using a semi-implicit 
Crank-Nicolson algorithm [381 ]: 



-1 _ j,n+2/3 _ 1 



h2/3 



o 2 S r 



h2/3 



1 



(17) 



in which the index i now labels the Fourier mode number. The derivatives <9 2 are discretized as in Eqn. [10] above. 
Note that Eqn. 1171 contains no mixing of the Cartesian components S Q/ 3 for any a[3 = xx, xy, yy ■ ■ ■ , so can be solved 



for each one separately. Bringing all terms in the unknown 23^ across to the left hand side, and putting all terms 
in the known 23™ +2 ^ 3 on the right hand side, we obtain a sparse set of linear equations characterized by a tridiagonal 



matrix on the left hand side. These are then solved for the 23^* using standard NAG routines [401 ]. 

Having updated 23" — > 23™+ 1 using the viscoelastic constitutive equation, we now update the stream- functions 
(f> n ,ip n — > <j> n+1 , ip n+1 using the x,y and z components of the force balance equation Again, we work in real 
flow-gradient space and reciprocal vorticity space. To eliminate the pressure from Eqn. 0] we subtract d y of the z 
component from d z of the y component to get the following equations, written separately for qi = and qi =/= 0: 



d 3 6\ n+1 



1 



dy^x 



n+1 

yl0(i+l/2) 



for qi = 0, 



(18) 
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0yC(£i/2) -^^C(ill/ 2 ) 



(19) 



f) i M n+1 — —f) V fnr « — fl 

a y^l0(j+l/2) — ?? a 2/ Zj 2/ z l0(i+l/2) 101 * — U ' 



(20) 



4„/,|n+l 



[i<7i3 y — £ M ) ™ +1 



-d 2 )£^ +1 ] for ft^O, 



(21) 



with V 2 = (dy—qf). The real and imaginary parts of these equations are treated separately. The third order equations 
(HHJ), (JTUJ) and (f2"0"|) are discretized at staggered half grid points 2/7+1/2 for j = 1 • • ■ A y — 1, with derivatives calculated 
as follows: 



dyf\j+i/2 = {fj+i - fj) 



(22) 



and 



(23) 



for any quantity /. (For clarity, the subscript i and the superscript n have been omitted from these expressions.) The 
fourth order equation is implemented at full grid points jjj, excluding those at the very edge of the domain (yo and 
t/jv). In it, dy and d v are discretized as in Eqns. [T01 and [Til respectively, and dy according to 



d 4 yf\ 3 



Ay' 



(/ 7+2 -4/ 7+1 + 6/ i -4/ j _ 1 + / j _ 2 ). 



(24) 



Each of Eqns. [18]to[2Tj for each mode index i, then takes the form of a sparse set of linear equations for , <>r 



These equations are solved using standard NAG routines [4C 



r 



It remains finally to specify the spatial discretization of 
the boundary conditions. In turn, this will prescribe the 
values of the flow variables on the phantom grid points 
that lie just outside the flow domain. 

In the vorticity direction z, the boundary conditions 
are periodic. For any quantity / on the grid zi • ■ • Zn z 
in real space Zi, we thus have f-i = In z ~i, fo = /w z , 
Jn z +i = fo, In z +2 = fi- In reciprocal space q h the 
periodic boundary conditions are always satisfied. 

In the flow-gradient direction y, the boundary condi- 
tions for the fluid velocity at the plates y = 0, 1 are those 
of no slip and no permeation. In terms of the stream- 
function </>, the no slip condition gives d y 4> = at y = 
and d y (j) = 7 at y = 1. We also note that (j> is only 
defined up to an arbitrary additive constant, and accord- 
ingly choose <j> = at y = 0. The third order Eqns. IT51 
and 1 191 then have three boundary conditions, as required. 
After discretizing real flow-gradient space yj in reciprocal 
vorticity space Qi, we then have 

<f>ii = 0, <t>io = 4>i2, and 0i(jv y +i) = 0i(JV B -i)+<S»o7Ay, 

in which Sij is the usual Kronecker delta function. 

In terms of the stream- function tp, the no-slip condi- 
tion gives d y ip — at y = 0,1, and the no-permeation 
condition gives <9 z V' = 0aty = 0,l. We also note that ip 
is only defined up to an arbitrary additive constant, and 



choose ij) = at y = 0. In the qi = equation (|2H)) . we 
then have 

tpn = 0, ipio = ipi2 and -0i(JV B +i) = ipi(N y -i)- 

These also hold for the qi ^ equation (|2ip , which obeys 
the additional condition 



iN, 



= 0. 



The zero-gradient boundary condition for the viscoelas- 
tic stress, after discretization on the flow-gradient grid 
2/1,2/2 -■•VNy, gives / i0 = fa and fi(N y +X) = fi(N y -i) for 
all components / = Ti xx , X™, ~S yy ■ ■ ■ . These apply in 
both real Zj and reciprocal % vorticity spaces. 

Extension to the case of a nonlinear grid in the y di- 
rection is straightforward in principle, but cumbersome 
in detail. Discretized derivatives are calculated via the 
usual Taylor expansions. For example, to first order accu- 
racy, centred second derivatives with respect to y become 



d 2 f- = 



Vj+i -2/3-1 



/j'+i - fj fj ~ fj- 



Vj+i ~ Vj Vj -Vj-i 



(25) 



We checked our nonlinear mapping carefully by perform- 
ing a few runs for identical parameter sets with both 
linear and highly nonlinear grids. All the numerical re- 
sults in this paper are converged with respect to grid and 
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FIG. 3: Steady state at a = 0.3, r) = 0.05, j = 2.0, £ = 
0.00375, L z — 4.0. Left: greyscale of T, xx m y — z plane. 
Middle: (y — z) velocity vectors, showing Taylor-like rolls. 
Right: vorticity banding of viscoelastic shear stress. 

timestep, to within the accuracy resolvable on the plots 
presented. 



VI. NONLINEAR STEADY STATE 

In each simulation run, we input as an initial condi- 
tion the ID gradient-banded state discussed in Sec. IIII1 
superposed with Fourier perturbations of tiny random 
amplitudes. As expected, under conditions where the lin- 
ear analysis of Sec. IIVI predicts the ID initial state to be 
unstable with respect to perturbations with wavevectors 
q = q z i in the vorticity direction, we find that these ini- 
tial disturbances grow in time. Full agreement between 
(i) the early-time growth rate and functional form of the 
fastest growing mode and (ii) the most unstable eigen- 
value and eigenfunction of the linear stability analysis 
provides a stringent check of our numerical method. 

During the instability, the initially flat interface be- 
tween the bands develops undulations that grow in time. 
At long times, once nonlinear effects become important, 
these saturate in a finite amplitude interfacial undulation 
to give a 2D steady state (Fig.[3]). The wavelength of the 
steady undulations corresponds to that of the maximally 
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FIG. 4: Selected stress: ID initial state and 2D steady state. 
L z = 2.0, a = 0.3, T] = 0.05, 7 = 2.0. At i = 0, the stress 
varies erratically about an average that depends on the grid. 
Inset: evolution to steady state at I — 0.00375. 

unstable mode of the linear analysis. For the parameters 
of Fig. [31 this is roughly twice the gap width. Associ- 
ated with these undulations are Taylor-like velocity rolls 
stacked in the vorticity direction (Fig. [3J middle), ac- 
companied by undulations of the stress along the wall 
(Fig. [3J right). The results of Fig. [3J could be tested ex- 
perimentally as follows. Optical measurements should re- 
veal birefringent stripes stacked in the vorticity direction, 
each of height comparable to the gap width. Likewise, 
the velocity rolls could potentially be measured using ve- 
locimetry. This is a challenging task, however, because 
the highest speed in Fig.[3fc is only 0(0.01). 

Our results capture recent experimental observations 
in which an initially ID gradient banded state of a worm- 
like surfactant solution was found to destabilise with re- 
spect to interfacial undulations with wavevector in the 
vorticity direction [I?]]. Indeed, several features of our 
results can be directly compared with these experiments, 
as follows. In the ultimate steady state, the wavelength 
O(L) and amplitude O(Z/10) of the undulations in our 
Fig. [3J are comparable to those in Fig. 2 of Ref. [27| . 
measured in units of the gap width L. The wavelength 
of these undulations was furthermore reported to increase 
with increasing average applied shear rate 7 [27| . consis- 
tent with the inset of our Fig. [2Jd. 

The kinetics of the instability can also be compared, 
via the temporal evolution of the stress signal. In the 
experiments [27|, a shear startup protocol was followed. 
Accordingly, the stress signal showed an initial overshoot 
followed by a decay (at 7 = 30s -1 ) on a timescale O(r) to 
a plateau value. This part of the dynamics corresponded 
to the initial formation of ID gradient bands. It is absent 
from our simulations, because we take as our initial con- 
dition an already ID gradient banded flow. Subsequently, 
the stress signal in Ref. [13] slowly increased by about 1% 
on a timescale O(100t). This part of the dynamics was 
associated with the ID gradient banded state destabilis- 
ing to exhibit vorticity undulations. As shown in Fig. [4l 
it is captured very well by our simulations: we find a slow 
stress increase 0(1%) on a time-scale O(100r), consistent 
with the experiments. 
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FIG. 5: Profile T, xx normal to the interface; £ = 0.00375 (thick 
lines), £ = 0.00250, 0.00125 (thin), a = 0.3, r\ = 0.05, j = 2.0, 
L z = 2.0. a,b) ID state with interfacial thickness d oc I. (c,d), 
(e,f) 2D state at the coordinate where the interface has 
maximal + and — y displacements. The thickness d appears 
independent of I, but note the bump of thickness 0(1). The 
offset yo (I) is chosen to centre the profile at the origin in each 
case. 

Some differences between our work and the experi- 
ments of Ref. [27j are noted as follows. In Ref. [27| . 
the instability was studied using light scattering tech- 
niques, which couple to concentration fluctuations. In 
the present manuscript, we do not consider concentra- 
tion coupling. In future work, it might be interesting 
to perform analogous simulations in the concentration 
coupled model of Ref. [4lj]. However, an important find- 
ing of the present work is that concentration coupling is 
not actually needed to trigger the basic undulatory in- 
stability. Indeed, we believe this to stem instead from 
normal stress effects, with concentration coupling a sub- 
dominant feature. Finally, we have not seen the exotic 
dynamics reported at the edges of the stress plateau in 
Ref. [27( . We cannot access small enough I for the insta- 
bility to persist here. (Recall Fig. [2Jd.) 



We return to comment on the validity of restricting 
our study to the y — z plane, which was done mainly for 
computational efficiency. As seen from our results, the 
ultimate amplitude of the interfacial undulation in this 
(y — z) plane is in fact comparable to that reported in 
the (x — y) plane in Ref. [34j (to within 10% at a = 0.3, 
T) = 0.05, 7 = 2.0, £ = 0.00375, L {x . z} = 2). The present 
study therefore shows that vorticity structuring is indeed 
important, but also that a full 3D simulation should be 
performed in future work. 

Finally, we discuss briefly the thickness d of the inter- 
face between the bands. In the ID initial state, d = 0(1). 
In the limit £ — ► 0, this gives an unphysically sharp in- 
terface d — > 0. Associated with this is a pathological 
steady state that strongly depends strongly on the flow 
history (42j. In 2D, in contrast, d appears virtually inde- 
pendent of £, as shown in Fig. [SJ>f. This is an important 
finding that could potentially obviate the gradient term 
P\7 2 S in Eqn.[3l which is needed to give a finite interfa- 
cial thickness in ID [l5| . Nonetheless, the interfacial pro- 
file does retain a small bump of thickness 0(1), Fig. [5^-f, 
suggesting that the local case I = remains pathological 
even in 2D. Indeed, at I — the stress signal varies errat- 
ically about an average that varies between runs, Fig. 31 
though purely numerical instability cannot be ruled out. 
This important issue will be pursued further in future 
work. 

To summarise, we have identified a mechanism by 
which vorticity stress bands and Taylor-like velocity rolls 
can form in a complex fluid, triggered by the instability 
of gradient shear banded flow with respect to interfacial 
undulations. In any real startup experiment, we would 
expect the vorticity instability to commence during the 
final stages of the initial band formation: above we as- 
sumed a complete separation of timescales between the 
processes. In future work, we will study the true dynam- 
ics of shear startup experiment in curved Couette flow. 
We will also extend to 3D, to study the interplay of vor- 
ticity banding with the dynamics of Ref. (34| . Robustness 
of the mechanism in other models will also be studied. 

The author thanks Peter Olmsted, Paul Callaghan, 
Georgina Wilkins and Helen Wilson for discussions, and 
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